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Abstract 

A recent claim that in quantum chromodynamics the gluon propagator vanishes 
in the infrared hmit, while the ghost propagator is more singular than a simple 
C^h! pole, is investigated analytically and numerically. This picture is shown to 

^ • be supported even at the level in which the vertices in the Dyson-Schwinger 

I equations are taken to be bare. The running coupling is shown to be uniquely 

• determined by the equations and to have a large finite infrared limit. 

■ 1 Introduction 

The proof of the renormalizabiUty of non-AbeUan gauge theories like QCD[||, and the discovery 
of ultraviolet asymptotic freedom |^], heralded a new phase in the acceptance of quantum field 
theories as serious candidates for the quantitative description of the weak, electromagnetic and 
strong interactions. Since the running coupling in QCD decreases logarithmically to zero as 
the renormalization point is taken to infinity, it seems reasonable to calculate it perturbatively 
in the deep ultraviolet regime, where it is very small, even though a proof is lacking that the 
perturbation series makes sense (for example, that it is strongly asymptotic). 

Although one is not sure that perturbation theory is reliable for QCD at very high energies, 
at very low energies it is quite clear that it is inadequate. Chiral symmetry breaking and 

fermion mass generation are typically non-perturbative phenomena. The obverse of ultraviolet 
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asymptotic freedom is infrared slavery or confinement. Since tlie coupling decreases as the 
energy increases, it increases as one goes to lower energies, and the possibility is open that its 
infrared limit is infinite. Many attempts Q — necessarily of a non-perturbative nature — have 
been made to show this divergence of the coupling in the infrared limit. Mandelstam initiated 
the study of the gluon Dyson- Schwinger equation in Landau gauge[0]. Although he did consider 
the gluon-ghost coupling, Mandelstam concluded provisionally that its effect could safely be 
neglected. This assumption was also made in subsequent workj^, A deficiency of these 
attempts to show that the gluon propagator is highly singular in the infrared is the necessity 
to posit certain cancellations of leading terms in the equations. An uncharitable case of petitio 
principii might almost be made (i.e. circularity). 

Recently, a new possibility has been opened up by the work of von Smekal, Hauck and Alkofer[^. 
In this work the coupling of the gluon to the ghost was not neglected. These authors claim that 
it is not the gluon, but rather the ghost propagator that is highly singular in the infrared limit. 
The running coupling itself has a finite though quite large value in the limit of zero energy, 
presumably large enough to guarantee chiral symmetry breaking in the quark equation [^]. 

In the present paper we investigate the claims made in the new work. We shall write the gluon 
propagator in Landau gauge as 

where a and b are colour indices, and where A = is the projection operator 

The ghost propagator will be written in the form 

G-\p) = -5'^^\g(-p'), 

and we shall refer to the scalar functions F and G as the gluon and ghost form factors, respec- 
tively. 

The claim made in Ref.Q is that, in the infrared limit x = —p^ — > 0, these form factors have 
the following behaviour: 

F(x)~x2^ (1) 
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where k ~ 0.92. To obtain these results certain Ansatze were made for the three-gluon and 
ghost-gluon vertices, functional forms inspired, but not uniquely determined by Slavnov- Taylor 
identities. In fact the Ansatz made for the ghost-gluon vertex is such that actually the infrared 
behaviour Eq. (||) is not consistent with the Dyson-Schwinger equations. The difficulty is the 
occurrence of a term 

'^^F{y)G\y) (2) 
Jx y 

in the equation for the ghost form factor, which, with the form (||), would yield an impermissible 
logj; factor in the limit x — > 0. Von Smekal et al. circumvent this problem by replacing one of 
the factors G{y) by G{x), thereby undermining to a large extent their supposed improvement of 
the vertex Ansatz. 

Since we found the ad hoc nature of this last replacement questionable, we decided first to 
see what would happen if one simply replaces the full vertices by bare ones. In this case the 
problematic logarithm of Eq. (||) does not occur, and we can simply analyze the equation as it 
stands. If the behaviour (||) were to go away, it would bode ill for the new approach. However, 
our finding is that, with bare vertices, the form (|l|) indeed remains good, but with the index 
changed to k ~ 0.77. Moreover, we can show that the solutions of the coupled gluon and ghost 
equations lie on a three-dimensional manifold, i.e. the general solution has three free parameters; 
nevertheless all solutions have the infrared behaviour (||). Our primary purpose in this initial 
paper is to explain the above findings in detail. 

In Landau gauge, the QCD Dyson-Schwinger equations lead to the following coupled integral 
equations for the renormalized gluon and ghost form factors: 

5' ^ Z-^' dq^ 



G-Hp')=Zs - r dqVGiq') r de'^Fir'), (3) 



dq 

P' 

—^Z,j^ dqqG{q)j^ dO-^ 
with r'^ = p^ + (f' — 2pq cos 9. The kernels are 

\ 2 p'' I 2 

Q[P,q,r) - \^^^+^P 4 +2+p2l^4 + lg2 2 2 +p2l^2 



Here the full three-gluon and the ghost-gluon vertices have been replaced by their bare values, 
while the four-gluon and quark-gluon vertices have been provisionally thrown away. To ob- 
tain these equations from the Dyson-Schwinger equations, we performed a Wick rotation and 
evaluated two trivial angular integrations. 

The form factors and the QCD coupling are renormalized using some renormalization prescrip- 
tion, Zi, Z3, Zi, Z3 being the renormalization constants for the triple gluon vertex, the gluon 
field, the gluon-ghost vertex and the ghost field defined by 

3/2 / — ~ 

F(/) = Z3F(p2) G(p2) = Z3G(/) 3 = ^50 = ^^50, (4) 

where Fip^), G{p'^) are the unrenormalized gluon and ghost form factors, -F(p^), G{p'^) the 
renormalized ones, is the bare coupling and g its renormalized value. One usually writes 
the renormalization constants and the renormalized quantities as functions of a renormalization 
scale /U, corresponding to a specific renormalization prescription. However, we will see in the 
following sections that the renormalization prescriptions can be made more general than is 
usually done in perturbation theory, and that each prescription will correspond to a solution of 
the non-perturbative integral equations. Instead of having the usual invariance of the running 
coupling with respect to a variation of the renormalization scale, we will find a more general 
invariance under an arbitrary transformation in the three-dimensional space of solutions of the 
integral equations. 

We wish to solve the coupled integral equations (|3|) for F and G, and we propose to do that in 
a future publication. For the moment we introduce a further simplification, the y-max approxi- 
mation. This amounts to replacing F(r^) and G{r^) in Eq. (^) by F{p^) and G{p'^) if > q^, 
but by F{q'^) and G{q'^) if < q'^. This approximation facilitates the analytical and numerical 
analysis of the equations, since the angular integrals can now be performed exactly, and indeed 
the resulting one-dimensional Volterra equations can be converted into nonlinear ordinary dif- 
ferential equations. This y-max approximation is very widely employed for these reasons, in 
particular by von Smekal et al., and although we do the same thing here, let us sound a note of 
warning: although we do not expect the qualitative picture of Eq. (||) to change, we do expect 
the value of the index k to be different when we treat the coupled equations without the y-max 
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approximation. We have already seen that k is sensitive to the choice of Ansatz for the vertex 
functions, and it is also affected by the y-max approximation. The bare vertex Ansatz is of 
course only a first guess; and it is clear also that the Ansatz of von Smekal et al. needs to be im- 
proved (the logarithm problem to which we alluded above persists when one no longer employs 
the y-max approximation). Nevertheless, the picture that von Smekal, Hauck and Alkofer have 
uncovered appears to be robust in its qualitative, and hopefully also in its semi-quantitative 
features: the gluon propagator is soft in the infrared (i.e. it vanishes in this limit, instead of 
blowing up like a pole), while the ghost propagator is hard (it is more singular than a pole). 
The consequences for the physics of the strong interaction need to be investigated. 

2 The coupled gluon-ghost equations 

The set of coupled integral equations for the gluon and ghost propagator, using the bare triple 
gluon vertex and the bare gluon-ghost vertex, and introducing the y-max approximation, is as 
follows: 



F-\x) = Z3 + XZ1 



+\Zi 



F{x) 



X \ x"^ 2x 
dy f 7y2 I7y 9 



X 2y 

A2 



(5) 



X 



2x2 2x 



F{y) + I 

Jx 



dy 
X y 



-7 + |)F^fe) 



Z3 — '^^^l 



F{x) '^y-G{y)+ '^F{y)G{y) 
'o X X Jx y 



(6) 



where A = g'^ /16tt'^, x = and y = ■ 

To solve Eqs. (^, we eliminate the renormalization constants Z3 and Z3 by subtracting the 
equations at 2; = o": 



F~\x) = F-\a) 

dy ( 7y2 I7y 9 



+\Z 



+\Zi 



F{x) 



X 



2x2 2x 



" dy 7y2 17y 9 



2cj2 2a 



dy n_x\ 
X y \8y, 



-7 r jFHy)+ I -[^.]F'\y) 



y \8yJ 



(7) 



F{y) 
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G-\x) = G-\a) 
--AZi 



(8) 



Fix) r'-^l Giy) - Fia) f dy) + T ^ 
X X Jo a a Jx y 



3 Symmetries of the reduced equations 



A very interesting simplification of Eqs. (|^, ^ is obtained if we throw away the gluon loop in 
Eq. (^), keeping only the ghost loop. This truncation is particularly interesting because, as 
we will show, its properties agree qualitatively with the requirements of a consistent physical 
picture, whereas the inclusion of the gluon loop introduces an ambiguity, which is probably due 
to the presence of terms involving the yet unknown renormalization constant Zi . The truncated 
set of equations is: 



F-\x) = F'^a) 



(9) 



G-Hx) 



G-\a)-^XZ, 



Fix) r Giy) -Fia) 

XX 



cr a Jx y 



(10) 



We will show that Eqs. (^, 10) have a three-dimensional space of solutions and that these 
solutions can be transformed into one another by means of simple scalings. 

First of all, if we have a solution Fix) and G(x), we can build a two-dimensional infinity of 
solutions simply by scaling these functions: 

Fix) = Fix)/a (11) 
Gix) = Gix)/b (12) 

which simply amounts to a redefinition of Z^ and Z^, i.e. to a change in the renormalization 
prescription. The new functions satisfy the same integral equations, with the rescaled coupling 
constant: 

A = Xab^ . 

Although the value of A is in general changed, this has no physical significance, since the following 
gauge invariant quantity is unchanged by the above transformations: 

AF(x)G2(x) = XFix)&ix). 
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Thus the two-dimensional manifold of solutions corresponds to the same physics, and one could, 
for instance, take A = 1 and set F{a) = 1 without essential loss of generality. 

A second, less trivial feature is the possibility to derive an infinite number of solutions starting 
from F{x) and G{x) just by scaling the momentum x to tx. The new functions F and G take 
the same values at momentum x as -F and G at momentum tx: 



Fix) 



F{tx) 



G{x) = G{tx) 



(13) 



In terms of the scaled quantities, 



we find 



X = x/t y = y/t a = a jt 



F-Ux] 



F-\a) 



G{x) ^ 

X 



3y 



2x 



3y 



^ + ^ Giy)-Gia) ^ -^ + -JL G{y)+ ^ G\y) 



cr 



2a- 



dy 



2y 



G-\x) = G-\a)-^\Zi 



Fix) r ^^G{y)-F{a) 

XX 



cr a Jx y 



This means that F{x) and G{x) are also solutions of the integral equations solved by F{x) and 
Again, all the solutions obtained by varying the scaling factor t correspond to the same 
physical picture, since a scaling of momentum merely corresponds to choosing the units for 
the momentum variable when renormalizing the coupling constant at a certain physical scale. 
It is clear that the three above-mentioned scaling properties allow us to construct the whole 
three-dimensional space of solutions starting from one specific solution. 

This three-fold scaling invariance has a physical relevance, as we will now show. From the 
renormalization of the gluon-ghost-ghost vertex, we define the renormalized coupling as 



a = Zi '^ZsZ^ ao , 



(14) 



where ao = 9o/^''^- According to Taylor 1^], Zi = 1 in the Landau gauge and, writing Eq. (14) 
with two different renormalization prescriptions, we find 
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aZ.^Z^' 



aZ:r^Z 



(15) 



Substituting Eq. (Q) for Z3 and Z3 and eliminating the unrenormalized quantities, we have 

aF{x)G^ix) = aF{x)G^{x) . (16) 



This ahows us to define the running couphng 

a{x) = aF{x)G'^ix), (17) 
which is independent of the renormahzation prescription. 

This invariance of a{x) with respect to the renormahzation prescription is exactly reproduced in 
the ghost-loop-only truncation. Furthermore, the renormahzation prescription is more general 
than what is normally used in perturbation theory. There the prescription usually is -F(/u) = 
G(/i) = 1 and a = a^"^^, and the invariance is taken to be an invariance with respect to the 
choice of the point /i at which these conditions are imposed. In our non-perturbative treatment 
the invariance is with respect to an arbitrary transformation in the three-dimensional space of 
solutions of the equations. For many of these solutions, there is even no such no such point at 
which F{^) = G{pL) = 1 or where a = ^ so that no traditional scale ^ can be attached 
to the renormahzation prescription itself: only the running coupling a{x) defined in Eq. ([T7|) is 
physically meaningful. We will see later that the loss of symmetry of the equations when we 
include the gluon loop in the current truncation scheme destroys this invariance with respect to 
an arbitrary renormahzation prescription, and different prescriptions no longer lead to the same 
physical running coupling. 

4 Infrared behaviour 

We will show analytically that the equations Eqs. (|^, ^) and Eqs. (^, |l^) have a consistent 
infrared asymptotic solution: 

F{x) = Ax^"" (18) 
G(x) = Bx-"", (19) 

and that these solutions even solve the ghost-loop-only equations (|^, |lO|) exactly for all momenta. 
Let us try the Ansatz 

F{x) = Ax"" G[x) = Bx^ . (20) 



8 



In the infrared asymptotic regime the gluon loop does not contribute to lowest order. Substi- 
tuting Eq. (^) into the integral equations |l^) we calculate 

"3 1 11 



22 + (3 3 + /3 4/3 



and 
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1 1 

+ 



2 + /3 a + P 



on condition that 



P> -2 



(21) 

(22) 
(23) 



to avoid infrared singularities. The powers on both sides of Eqs. (21, 22) agree if 



and defining the index k by 



a = -2/3, 



a = 2k j3 = —K 



we find that both the constant and the power terms in Eq. (21) and Eq. ( |22D match if 



XZiAB^ 



and XZiAB'^ 



1 1 

+ 



T -1 



2(2 - k) 3 - k 4k 



1 1 



2 — K K 



(24) 



(25) 



(26) 



Elimination of XZiAB^ yields a quadratic equation for k, which remarkably does not depend 
on the value of the coupling strength A: 



19k^ + 77k + 48 = . 



which has two real solutions 



77 ± V2281 
38 ' 



or 



Ki PS 0.769479 and Kg ^ 3.28315. 



(27) 

(28) 
(29) 



The second root is spurious: it must be rejected because it gives rise to infrared singularities 
and thus does not give a solution of the integral equation. 

Replacement of k, by ki in Eq. ([25|) or Eq. (p^) yields the condition: 



V = XZiAB^ w 0.912771 . 



(30) 



From Eq. ([T^) we know that the running couphng is given by 



a{x) = 4ttXF{x)G^ 




(31) 



in the Landau gauge. Condition Eq. (^) is important, as it tells us that the running coupling 
has a non-trivial infrared fixed point 



This means that the ghost field, which only introduces quantitative corrections to the pertur- 
bative ultraviolet behaviour of the running coupling, does alter its infrared behaviour in a very 
drastic way. 

We will show further on that the running coupling remains almost constant up to a certain 
momentum scale x, after which it decreases as 1/logx. The momentum scale at which the 
constant bends over into a logarithmic tail is closely related to the value of Aqcd- This is easily 
understood intuitively, since the perturbative ultraviolet behaviour of the running coupling blows 
up very quickly as the momentum gets down to 0{Aqcd)- 

5 Infrared asymptotic solution 

Although we have seen in the previous section that the pure power behaviours for F{x) and G{x) 
solve the reduced equations exactly, these power solutions only give rise to a two-dimensional 
space of solutions. However, the numerical results told us that the equations were much richer 
then we initially believed. These numerical results tended to suggest that the power solutions 
are only one very special two-dimensional family of solutions in the midst of a whole three- 
dimensional space. Typical non-power solutions showed an infrared behaviour completely con- 
sistent with the power solution mentioned earlier, which then bends over quite rapidly at some 
momentum x into a completely different ultraviolet behaviour which seemed to be proportional 
to some power of the logarithm of momentum. A straightforward investigation of the ultraviolet 
asymptotic behaviour of the solutions tells us that such powers of logarithms are indeed con- 
sistent ultraviolet solutions, but no obvious mechanism seemed available to match the infrared 
to the ultraviolet parts of the solutions, making us believe at first that the numerical program 
was giving us spurious pseudo-solutions, due to some numerical inaccuracies or artifacts. One 
of the main reasons was that the infrared power behaviour only contains one free parameter, 



lim a(x 

x-*0 




11.4702. 



(32) 
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and a standard asymptotic expansion does not add any corrections to the leading power. If the 
infrared asymptotic solution contains only one parameter, it was very unclear how an infinite 
number of solutions with log-tails could develop out of each power solution. Nevertheless the 
numerical results indicated that each power solution had an infinite number of corresponding 
log-tailed solutions, and each solution seemed to be characterized by the momentum at which 
the log-tail sets in. 

The traditional asymptotic expansion one would normally try, is as follows: 

N 

F{x) = x^^^Aix' (33) 

N 

G{x) = x-^'^BixK (34) 

i=0 

The reason for this is that each term in the expansion usually generates terms, through integra- 
tion, that are of the same power or one unit higher. However, the fact that the equations under 
consideration are exactly solved by the power solution alters the reasoning. The leading power 
term does not generate additional, next-to-leading order terms, and all Ai,Bi for i > have to 
be zero for consistency reasons. 

However, the fact that the power solution solves the integral equations does not mean that this 
is the unique solution, and we next tried an infrared asymptotic solution of the shape: 

F{x) = Aqx^"" + Aix""' (35) 
G{x) = Box^'^ + Bix^K 

with ai > 2k and /3i > —k. Substitution of these solutions into Eqs. (P, ^), tells us that 
consistency is obtained if ai — Pi = k, as for the leading power, but it gives an additional 
constraint, fixing the value of the exponent of the next-to-leading exponent. However, the 
solution proposed above does generate additional higher order terms, and consistent asymptotic 
infrared expansions can be built as follows: 

AT 

F{x) = x^^Y.^i^''' (36) 

i=0 
N 

G{x) = x-^Y.Bix'", 

where the exponents of successive powers always increase by the same amount p > 0. To 
check the consistency of these infrared asymptotic expansions, we substitute them into Eqs. (P, 
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10). We make a Taylor expansion of the left-hand sides of these equations and expand the 
series multiplications, before integration, on the right-hand sides. Consistency requires that the 
coefficients of equal powers of momentum match each other on both sides of the equations. 

The conditions on the leading term remain unchanged as described in Sect. ^, with k ~ 0.769479 
and u = XZiAqBq ~ 0.912771. Equating the second order terms on left and right-hand sides of 
both integral equations yields the following set of two homogeneous linear algebraic equations 
for ai = Ai/Aq and bi = Bi/Bq: 

ai / 3 1 3 1 1 \ , 

— + r h r ]bi=0 

V \2{2-K^p) 3-K-Fp 2(2 -k) 3-k -2k^p) 

1 1 \ /I 1 4 \ , 

ai + — \ 7. — + 7^ M^i = 



K -|- /9 2 — kJ \k + p 2 — k + p 9u 

This set of equations will only have non-trivial solutions if its determinant is zero, in which case 
it will have a one-parameter infinite number of solutions. The characteristic equation is: 

- 9.27685 p"^ - 15.5544 p^ + 30.2899 + 71.5686 p = 0. (37) 

The four solutions are: 

p = 0, p= 1.96964, = -1.82316 ±0.770012 i. (38) 

The solution p = corresponds to the pure power solution. The two complex solutions are 
spurious as they are not consistent with Rep > 0, while the solution p = 1.96964 gives rise to 
consistent infrared asymptotic expansions. 

The linear homogeneous set of equations then yields 

= 6i/ai = 0.829602, (39) 

and the solutions of this set of equations can, for instance, be parametrized by ai. 
Let us define 

an = An/Ao , bn = B^/Bq , (40) 

in terms of which we find the following heterogeneous set of equations for 02 and 62: 



V 



3 13 11 

+ 



2(2 -k) 3-k 2(2-fi; + 2p) 3 - k + 2/9 -2k + 2p, 



(41) 
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a 



1 



1 



1 



K + 2p 2 — K 



a2 + 



u \2{2- K + p) 3-K + /3 2(-2K + 2p) 
1 1 4 



hi 



+ 



K + 2p 2 — K + 2/9 yz^ 
_ 46f / 1 1 



&2 



K + 2p 2 — K + p 



with unique solution 



and for 03, ^3: 



02 = 0.408732 af 



62 = 1.31169 af 



as 



+ 



1 



+ 



1 



1 



2(2 -k) 3-k 2(2 -k + 3/9) 3-K + 3p -2K + 3p 



63 



2aia2 — af 



1 



1 



1 



1 



1 



K + 3p 2 — K 

^ 4(2bib2 - bj) 
9u 



2{2-K + 2p) 3- K + 2p2{2- K + p) 3- K + p -2k + 3p 
1 14 

as + 



(42) 
6162 



K + 3p 2 — K + 3p 9i/ 
1 1 



K + 3p 2 — K, + 2p 



ai&2 



1 



1 



K + 3p 2 — K + p 



0261 . 



with unique solution 



03 = -0.761655 af 



63 = 0.783905 a? . 



By induction one can prove that the higher order terms all yield sets of equations of the same 
nature as Eq. (^) and Eq. (^), where the right-hand side of the set defining the coefficients 
an,bn are proportional to a". This means that we have a general solution for the nth order 
coefficient of the type 

On = /nfll bn = ^nOi (43) 

for n > 1, where the fn,gn are constants (independent of A and of Zi). 
The asymptotic expansions Eq. (|3^ ) can thus be written in the form 



F{x) 
G{x) 



(44) 



where Bq and ai = Ai/Aq are chosen to be the free parameters spanning the whole three- 
dimensional space of solutions of Eqs. (^ [lO| ) in the infrared region, and where (to 6 significant 
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figures) 



V = XZxAqBI = 0.912771 , k = 0.769479 , p = 1.96964 (45) 

/i = l, /2 = 0.408732, /3 = -0.761655, 

51 = = 6i/ai = 0.829602, 52 = 1-31169, c/3 = 0.783905 , 

It is precisely the existence of a third independent parameter, namely oi, which allows the 
infrared power solution to bend over in a logarithmic tail in a way consistent with the integral 
equations. To build a solution that is both consistent with the infrared asymptotic expansion 
set up in this section and the asymptotic ultraviolet logarithmic behaviour which will be derived 
in the next section, the parameter oi has to be negative, as has been inferred from the numerical 
results calculated with the Runge-Kutta method and with the direct integral equation method. 
If ai = we retrieve the pure power solution and if ai > there does not seem to be a 
singularity- free solution for x E [0,A^]. 

As we have shown in Sect. |3|, the three-dimensional family of solutions can also be constructed 
once we have found one solution, just by relying on the three distinct scale invariances (^, |l^. 



13). How these scale invariances correspond to choices of infrared asymptotic parameters will 



now be elucidated. 



The function scalings ( |l l| , 12) of F{x), G{x) correspond to similar scalings of Aq,Bq in the 
infrared expansions Eq. (^4|), 



^0 = A^/a , Bq = Bo/b , 

such that condition Eq. (45) remains satisfied with A = Aa6^, and oi is left unchanged. 
Less trivial is the momentum scaling invariance of the space of solutions: 

F{x) = F{tx) G{x) = G{tx) . (46) 



Using these definitions in Eq. ([44D, we find, after some rearrangement. 



F{x) = {t'-A^)x^-l^ + Y,J,{tPairx'P^ 
G{x) = {t-''Bo)x-''l^ + Y,^gi{tPairx'P^ 
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This shows that the infrared expansions for the momentum scaled functions F{x),G{x) corre- 
spond to asymptotic expansions parametrized by 



Bq =t '^Bq 



and 



(47) 



and that the asymptotic expansions indeed obey Eq. (44) and the conditions Eq. (|4q). As 
we expected, there is a one-to-one correspondence between the solutions constructed from the 
scaling invariances based on the symmetries of the equations, and the parameters Aq, Bq and 
ai characterizing their infrared expansions. 

Let us now construct the asymptotic expansion of the running coupling (with Zi = 1) using the 
expansions (^^: 

/ N \ 



X{x) = XF{x)G\x) = XAoBI 



1 + 



iCLlX 



i=l 



N 



1=1 



or (again truncating at N) 



where 



A(x) = V 



(48) 



(49) 



hi = 2.65920 , h2 = 5.37956 , = 6.97232 , 

which tells us that the running coupling only depends on the dimensionful parameter ai = 
Ai/Aq^ and is independent of A, and Bq. Furthermore, we can show from Eqs. (^, ^) that 
the running coupling corresponding to the parameter oi, is identical to the running coupling 
with parameter ai after scaling the momentum with a factor t = (ai/aiY^^ . This tells us that 
the momentum units of ai are unambiguously related to the physical scale of the experimentally 
determined running coupling. 

We now introduce a momentum scale 

1 



{hi\ai\f'' ' 



(recall that ai < 0), such that 



N 



i=l 



(50) 



(51) 



where we defined 

hi 



hi 



h\ 



hi = 1 



h2 = 0.760753 , 



/i3 = 0.370785 , 



15 



We will see from the numerical results that is a good estimate of the scale up to which the 
infrared asymptotic expansion remains valid. 



6 Ultraviolet behaviour 



We now turn to the investigation of the ultraviolet asymptotic behaviour of the solutions. As 
discussed before, the numerical results show a three-dimensional space of solutions, which has 
been confirmed by an analytical study of the global symmetries of the integral equations and 
by the study of the infrared asymptotic expansions of the solutions. Except for the pure power 
solution, all these solutions bend over in a log-tail above a certain momentum scale x. We will 
now check the consistency of such ultraviolet logarithmic solutions. 

Suppose the solutions for F{x) and G{x), taking on the values and at some momentum 
^ in the perturbative regime, have the following ultraviolet behaviour: 
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F{x) 
G{x) 



F,, 



log ( - ) + 1 



UJ log 








+ 1 









(52) 
(53) 



We check the consistency of these ultraviolet solutions by substituting these expressions in 
Eqs. d, 0). 



The ghost equation ( [ToD yields, to leading log, 
-5 



g: 



wlogi - ) +1 



G- 



a 



L^log(-) +1 



After evaluating the integral we get 

-5 



UJ los 



+ 1 

4w(7-|-(5 + 1) 



UJ loe 



UJ log 



9 
4' 



r dy 


UJ log 


m 


L y 







a 
A* 

+ 1 



+ 1 

7+5+1 



7+5 

(54) 
(55) 



a; log ( - ) + 1 



7+5+1" 



Matching the index of the leading powers of logarithms in Eq. ( ^5| ) one finds the consistency 
condition: 

7 + 25 = -l (56) 
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and, equating the leading log coefficients in Eq. (55), using Eq. (|56[), we get 



(57) 



Substituting the solutions Eqs. (|52| , |53|) in the gluon equation Eq. and keeping only the 
leading log terms, we find 



wlogi -1+1 



cjlogi - I +1 



a;log( ^ ) +1 



25 



■ (58) 



After performing the integrals, we find 



log ( - ) + 1 



+ 



-7 



2a;(2(^+ 1) 



cjlog ( -j +1 

\ -1 25+1 

cjlog ' j^j + 1 



(59) 



wlog( -1+1 



l2<5+l" 



Consistency of the exponents on both sides of the equation is automatically guaranteed by 
Eq. ([5^). Then, equating the coefficients of the leading log contributions of Eq. (59), and 
substituting Eq. (p§), we obtain 



(60) 



From Eqs. (56, 57, 60) we then find 



1 

7 = o 
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and the equivalent conditions Eqs. (E% 30) yield 



(61) 



(62) 



Thus, the ultraviolet solutions for F{x),G{x) can be written as 



Fix) 



F„, 



4XZiF^Gilog{-] +1 



Gix) = G^ 



AXZiF^Gpogl-] +1 



-9/16 



and the renormalization group invariant running coupling is given by 

A(x) = XFix)G^{x) ^ 



(63) 
(64) 

(65) 
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We can rewrite this in the form 



A(x) = y ^ , (66) 

/3olog 







where /3o = 4, and the QCD scale is given by 



A^Ci. = /^exp(^-^3^j , (67) 

if Z\ = 1. We see that fixing XF^G^^ at a scale n, in the perturbative regime, indeed amounts 
to a definition the value of Aqcd- 

The leading-log coefficient is /?o = 4, but this is not in agreement with perturbation theory, 
where Pq = 11. However, the reason for this is obvious, as we only considered the ghost loop 
and discarded the gluon loop in the gluon equation. 

7 Numerical method 

Knowing the infrared and ultraviolet asymptotic behaviours of the coupled equations Eqs. (0, 



10), we now go on to solve the equations numerically in order to see if we can find consistent 
solutions over the whole momentum range, connecting both asymptotic regions, hopefully giv- 
ing us more insight into the transition from the regime of asymptotic freedom to the state of 
confinement. We first give a short overview of the numerical method we used. 

We use a numerical method developed by one of us for the study of dynamical fermion mass 



generation in QED4[10|. This method directly solves the coupled integral equations by an iter- 
ative numerical scheme. We also checked the results, found with the integral equation method, 
with a Runge-Kutta method applied to the set of differential equations derived from the integral 
equations. 

We now give an outline of the main features of the integral equation method. Unlike most other 
methods used thus far, we replaced the widely used discretization of the unknown functions 
by smooth polynomial approximations, introducing Chebyshev expansions for the gluon and 
ghost form factors F{x) and G{x) and using the logarithm of momentum squared as variable. 
To improve the accuracy of the Chebyshev approximations we first extract the infrared power 
behaviours of the form factors, although this only has a minor influence. The form factors are 
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approximated by 



F(x) 



Ax 



2k 



Gix) = Bx-^ 



N-l 

h 



with 



s\x] = 



logio(a;/Ae) 



(68) 
(69) 

logio(A/6) ' 

and where A is the ultraviolet cutoff, and e is the infrared cutoff, only needed for numerical 
purposes. We require both integral equations to be satisfied at N fixed external momenta, 
in order to determine the 2 Chebyshev coefficients Oj , 6j . Using smooth expansions has the 
advantage of allowing us an absolute freedom in the choice of quadrature rules used to compute 
the various integrals numerically. This is required if we want to achieve a high accuracy in our 
results. The integration region is first split into an analytical integral over [0, e'^] and a numerical 
integral over [e^,A^]. The integral over [0, e^] is computed analytically from the asymptotic 
infrared behaviour discussed in Sect. ^. This is needed as the infrared part of the integral 
is highly non-negligible, especially in the case of the gluon equation. For an efficient choice 
of quadrature rule we split the numerical integral into three regions, these are [e^, min(2;, cr)], 
[min(x, cj), max(x, cj)] and [max(x, cr), A^], where x is the external momentum and a is the 
subtraction point. The splitting of the region of numerical integration into three subregions 
is needed as the integrands are not smooth at the boundaries of these regions and too much 
accuracy is lost if one uses quadrature rules spanning these boundaries. A sensible choice 
of quadrature rule on each integration region is for instance a composite 4-points Gaussian 
integration rule, where the composite rules are delimited by the region boundaries and the values 
of the external momenta at which we require the integral equation to be satisfied. This setup will 
yield 2N coupled, non- linear, algebraic equations for the 2N Chebyshev coefficients Oj and hy 
In traditional Dyson-Schwinger studies, the unknowns are usually determined by what is often 
called the natural iteration method, where the current approximation to the unknowns is used 
in the integrals of the right-hand side of the equations in order to provide a new approximation 
to the unknowns used in the left-hand side of the equations. This iteration method however is 
not necessarily convergent, and when it is convergent it often converges very slowly, as has been 



shown in Ref.[10|. This slow rate of convergence is not only inefficient, but more importantly 
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it makes it very difficult to get a reliable estimate of the accuracy of the solution. For this 
reason, our numerical method uses the Newton method to solve sets of non-linear equations. 
This method uses the derivatives of the equations with respect to the unknowns to speed up the 
convergence. If the starting guess to the unknown coefficients is close enough to the solution, the 
convergence rate is even quadratic. Let us symbolically rewrite the coupled functional equations 
as follows: 

f{x)[F,G] = 
gix)[F,G] = 0, 

where / and g are equivalent to the Eqs. (^, |10|) and x £ [0, A^]. For the numerical solution, we 
require this equation to be satisfied at the external momenta Xi, the functions F and G are ex- 
panded as Chebyshev polynomials with coefficients Oj and bj, and the integrals are approximated 
by a suitable quadrature rule. The equations then become 

f{xi)[aj,bj] = 
g{xi)[aj,bj] = 0, 

where i,j = . . . N — 1 and / and g are the numerical approximations to / and g when the 
integrals are replaced with quadrature rules. 

The Newton method will yield successive approximations to the solutions, given by 

and the {n + l)-th improvements Aa""*'^, Ab^^^ are given by the solutions of the INxlN set of 
linear equations 

ottj ■' obj ■' 

oaj ■' obj ■' 

where the equations are taken at the external momenta Xi and each equation includes implicit 
summations over j. 

The total accuracy depends on the combination of the accuracies of the Chebyshev expansion 
and of the quadrature rule and on the convergence criterion of the Newton iteration. We will 
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see later, after comparing the results of this method with those obtained with the Runge-Kutta 
method, that we can achieve a very high accuracy over quite a broad momentum range. 



Using this method, we performed a meticulous study of the equations Eqs. 10). We note 
that, for a fixed value of A, the equations have two free parameters, for instance F{a) and G{a) 
(restricted by Eq. (|30|) , XF{a)G'^{a) < 0.912771). Furthermore, as shown in Sect. a scaling 
of A can always be absorbed in a redefinition of the unknown functions F{x) and G{x), such 
that knowing the solution space for one value of A, we can build the solutions for any arbitrary 
value of A in a straightforward way. Moreover, such scalings of A leave the running coupling 
XF{x)G'^{x) unchanged. 

In practice, we choose an alternative pair of parameters, F{a) and A, where A is the leading 



infrared gluon coefficient defined in Eq. (18). The choice of these two parameters is suggested 
by the numerical solution method. Using Eq. (|30|), the value of A also determines the leading 
infrared ghost coefficient B, and allows us to compute a quite accurate analytical approximation 
to the infrared part of the integral over [0, e^], if is sufficiently small. The choice of F{a) as 
second parameter can be viewed as a measure of the deviation from the pure power behaviour 
at momentum a. We have taken the subtraction scale to be a = 1 and varied both parameters 
A and F(l) = Fi for a fixed value of A = 1 and Zi = 1. 



8 Results 



We vary both parameters A and Fi to scan the two parameter space of solutions, keeping A = 1 
fixed. As expected we retrieve the scaling invariances discussed in the previous sections. 

If we plot the solutions for F(x) and G{x) for various sets {A, Fi), as in Fig. |^, we can check 
that every solution can be transformed into another one by a unique transformation (t,r) cor- 
responding to a momentum scaling tx and a function scaling rF{x), G{x)/y/r. The numerical 
results clearly show the expected power behaviour in the infrared region and the logarithmic 
behaviour in the ultraviolet region. The value of the exponents and of the coefficients in front of 
the power in the infrared region is completely consistent with the analytical treatment of Sect. ^, 
which was also used to compute the infrared part [0, e^] of the integrals analytically. As can be 
seen from the plots, the gluon form factor, which starts off as a power with a given coefficient 
A, will bend over at some cross-over point x, such that the further logarithmic behaviour of the 
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Figure 2: Running coupling a{x) = aF{x)G'^{x) versus momentum x (on log-log plot), for A = 1 
and Fi = W~^{a), lO-'^(b), lO'^^c), 0.01(d) and 0.1(e). 



function consistently leads to a value Fi at the subtraction scale a = 1. The logarithmic be- 
haviour of F{x) and G{x) also satisfies the ultraviolet leading log behaviour analyzed in Sect. 
It is remarkable that both asymptotic regimes, infrared and ultraviolet, seem to connect onto 
each other at some momentum x, with scarcely any intermediate regime. 

If we look at the running coupling we see that all the solutions are just translations of each other 
when plotted on a logarithmic momentum scale, as is illustrated in Fig. |2|. This corresponds to 
the invariance of the space of solutions with respect to scaling of momentum. It also shows the 
physical equivalence of all solutions as such a transformation can always been absorbed into a 
redefinition of momentum units. 

It is also interesting to compare the numerical results with the analytic asymptotic calculations in 
order to investigate in which momentum regions the asymptotic solutions are valid. As example 
we consider the case A = 1 and Fi = 0.1. To compute the infrared asymptotic expansion we 
need to know the value of the infrared parameter oi in Eq. (^^. We used the Runge-Kutta 



method, which will be described in Sect. IC, in order to determine the value of ai yielding a 
value of Fi{l) = 0.1 for the gluon form factor, with A = I. For this specific case, the value is 
ai ~ —10.27685 or Q"^ ~ 0.186475 (from Eq. (|50[)). The infrared asymptotic expansion, derived 
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Figure 3: Running coupling aF{x)G'^{x) versus momentum x (on log-log plot), for parameter 
values A = 1 and Fi = 0.1 together with its infrared asymptotic expansion and its ultraviolet 
asymptotic behaviour. 



in Sect. |5[ is calculated from Eq. (^) and truncated after four terms: 



a{x) ~ 47rz/ 



iL-)\ 0.760753 (i)"- 0^370785 (i)" 



The ultraviolet asymptotic behaviour, derived in Sect. |6|, is described by Eq. (| 



a{x) 



4 log 







where we use Eq. 



^QCD = /^exp 



to compute the value of Aqcd for the case under consideration. We choose fi in the perturbative 
regime, for example /U = 1032.15, where the numerical results yield A(/i) = AF(^)G^(/i) ~ 
0.0259676, and find 

A|cD = 0.06802, 

still in arbitrary units. 

In Fig. |3| we plot the running coupling versus momentum, together with its infrared asymptotic 
expansion and the ultraviolet asymptotic behaviour. The agreement between the analytical and 
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numerical results is extremely good, and it can be seen that both asymptotic behaviours flow 
into each other, almost without any intermediate regime. The vertical line in Fig. ^ situates the 
scale of Aq(^£). We see that ^^CD ^^^^ momentum regime where the infrared asymptotic 

expansion has already taken over from the logarithmic behaviour, and where the running cou- 
pling has become almost constant. Furthermore, the infrared scale Vt^ ~ 0.186 seems to be a 
good measure to delimit the infrared region where the asymptotic expansion is valid. 

We can even give a numerical relation between and ^^CD (where the latter is computed 
from leading log only), namely 

and the ratio is the same for all solutions of the equations ^). The simple relation between 
Vt and Aqcd is a consequence of the symmetries of the ghost-loop-only truncation. If we include 
the gluon loop in the same truncation scheme, the asymptotic expansion of the running coupling 
will no longer depend on oi alone, but on the other parameters of the infrared expansions as 
well. Hence an ultraviolet renormalization, leading to a specific value of Aq^^, will correspond 
to a family of running couplings, all having slightly different behaviours in the intermediate 
regime, and there will be an ambiguity in the determination of the non-perturbative running 
coupling. 

The units in which Aqcd is expressed are still arbitrary, as we still have to compare the nu- 
merical results with experimental data. In the truncation under consideration, comparison with 
experimental results is only interesting to check our methodology. The numbers which will 
come out of the analysis are not to be taken too seriously, as the leading-log /3-coefficient in the 
ghost-loop- only case is 4, while in the pure gauge theory it ought to be 11, and in the presence 
of flavours of fermions it should become (33 — 2nj)/3, as is known from perturbation theory. 
This means that the ultraviolet running of the coupling will be too slow in the case we are 
considering, and the value of Aqcd will come out far too low. If we use the PDG|jll| value 
as{Mz = 91.187 GeV) = 0.118, we find A^^.^ = 2.277 x 10"^ GeV^ (using Eq. Thus, 
if we want to map the numerical results discussed above (expressed in arbitrary momentum 
units amu) to physical reality we have to set l{amu)'^ = 3.348 x 10^^ GeV^ and the mass of 
the Z-boson will be at 2.484 x 10^'^ (amu)^. Of course all this is only hypothetical as the (3- 
coefficient of the physical theory should be 25/3 in the presence of 4 fermion flavours, and this 
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would yield a value t^Qjy = 0.02343 GeV^. We also note that the freedom of renormalization 
allows us to choose the renormalized quantities, as is usually done in perturbation theory, i.e. 
as{Mz = 91.187 GeV) = 0.118, F{Mz) = G{Mz) = 1. This is possible, as all solutions of 
the three-dimensional space of solutions are physically equivalent. We will see later that this 
is not as obvious as it seems, and that it depends on the truncation scheme. Including the 
gluon-loop (with Zi = 1) in the gluon equation, either with a bare triple gluon vertex or even 
with a Ball-Chiu vertex, will destroy this invariance and different solutions will correspond to 
couplings running in different ways. 

9 Starting guess 

The Newton method, which is at the core of our numerical method, is a quadratically convergent 
iterative method, if the initial approximations to the unknown functions are sufficiently close 
to the exact solutions. The meaning of sufficiently close depends however entirely on the kernel 
of the integral equation. We observed that for the coupled gluon-ghost equations, the starting 
guess must not be too remote from the exact solution, if the method is to converge. This is 
in contrast with previous work on chiral symmetry breaking in QED and on the Mandelstam 
approximation to the gluon propagator in QCD, where the method was extremely insensitive to 
the starting guess. 

It turns out that in the case of the coupled gluon-ghost equations, the starting guesses have to 
be chosen quite sensibly, especially in the asymptotic regions. In practice, we used the analytic 
asymptotic solutions to build good enough starting guesses for the form factors. 

Given the parameters A and Fi, the leading order infrared ghost coefficient is 



which can be seen as a crude approximation to the bend-over point. A possible construction is: 




and we define x as 




-1 2k ^ 
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-1 -9/16 



G{x) = B 
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Figure 4: Comparison of the solutions for F{x) and G{x) with their starting guesses used in the 
iterative Newton method, for A = 1, ^ = 1 and Fi = 0.1. 

which has the correct leading infrared asymptotic behaviour for F{x) and G{x) and agrees well 
with their leading ultraviolet logarithmic behaviour, as is illustrated in Fig. ^. 

Although it seems that the starting guesses F{x) and G{x) are extremely close to the eventual 
numerical solutions, we see that the Newton method does alter the running coupling a(x) = 
47rAF(x)G^(x) substantially while converging to the solution, as is shown in Fig. |5[ 

10 Runge-Kutta method 



Rewrite the equations (|^, |T^), for cr = 1, as 

-G{x 



F-i(x) =r/ + A 



x-^ Jo \ 2 X Jx 2y 



and 



G-Hx) = c-h 



Fix) 



Jo 



r dyyG{y)+ j' '^F{y)G{y) 
Jo Jx y 



(71) 



(72) 



where rj and C are constants. As discussed before we can choose A = 1, since an arbitrary value 
of A can be recovered by applying an appropriate scaling to the form factors F{x) and G{x). 

Let us rewrite the above equations in the form 

3^_.._ 1 f^dy^2^ 



F-\x) = rj + -G{x)K{x) - G{x)L{x) + -f -^G\y) 
^ ^ Jx y 



(73) 
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Figure 5: Comparison of the solution for the running couphng a{x) with its starting guess, for 
A = 1, ^ = 1 and Fi = 0.1. 



and 



where 



and 



G-\x) = C - If{x)K{x) ^F{y)Giy) 

K{x) = ^ rdyyG{y) 
x'^ Jo 



L{x) = ^ rdyy^G{y). 

X JQ 

On differentiating the above four equations, we obtain 



F 
G 



F^l-^GK - \Gk + GL + GL\ + iF^G^ 



IG'[FK + FK] - IFG'' 
k = G-2K 



(74) 
(75) 
(76) 



L = G-3L, 

where F = ^ = xj^ etc., with t = logx. After a Httle algebra, we can throw the first two of 
these equations into the form 



F = 3F{X -Y)-FZ{^X -Y) 
G = ZG, 



(77) 
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where 



X = FGK 



Y = FGL 



X(3X - 3y - 2) 

^ - | + x(|x-y)- ^'^^ 

This system is suitable for an apphcation of the Runge-Kutta method; but we must first address 
the question of the existence and multiphcity of the solutions, first of the differential system, 
and then of the integral equations ( |73| , ^) from which they were derived. We already know an 
(exact) solution, namely 

F{x) = Ax^"" 

G{x) = Bx-^, (79) 
where k was defined in Eq. (p^), and where (compare Eq. ( p6|) with A = 1) 
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A priori we would expect there to be four free parameters for the differential system, corre- 
sponding say to the values of-F(l), G(l), K[\) and L[l)^ from which the differential equations 
could step-by-step be integrated, for example by the Runge-Kutta method. In general, solutions 
of the differential system would not satisfy the requirements x^K^x) and x^L{x) ^ as 



X — > 0. In fact the lower limits of the integrals in Eqs.(^, 76) would be incorrectly replaced by 
nonzero constants. Imposing the requisite boundary conditions at x = 0, we expect to reduce 
the number of arbitrary constants in the general solution from four to two. Since there is a scal- 
ing invariance that leaves FG^ unchanged, this means that, after we have removed this trivial 



degree of freedom by fixing A in Eq. (79), we should still have one non-trivial free parameter. 
Where is it? 

In Sect. P we have seen that we can construct the following infrared asymptotic expansion, 
Eq. (^^, for the general solutions F{x) and G{x): 



F{x) = Aox^^l^ + j^Ua'i^''^ 
-- Box-^l^ + Y^^QiO^^x"^ 



(81) 



N 

G{x 
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We expect these series to have zero radius of convergence, so they have been truncated in the 
anticipation that they are asymptotic series — that is, for smaU values of x, there wih be an 
optimal truncation point, N, for which the finite series is a good approximation. From this 
expansion we see that, besides the parameters and Bq, there is one more free parameter, ai. 




On substituting the series for G in the definitions Eqs. (75, 76), we find 

Kix) = Box~^(-^ + Y: ^ (82) 

L{x) = 

The knowledge of the infrared asymptotic expansions for F{x), G{x), K{x) and L{x) allows us 
to use a Runge-Kutta method, starting from a momentum point deep in the infrared region 
and building the solution for increasing momenta. The Runge-Kutta method was run using the 
NDSolve routine of Mathematica 3.0. The problem is solved as a function of t = logx and as 
the starting point, the IR series Eqs. (81, [8^ ) are evaluated at x = 0.0001 with N = 8, using the 
coefficients fj and gj which are calculated with Mathematica as well. The Runge-Kutta routine 
is run with 25 digit precision and 10,000 steps from x = 10~^ to x = 10^ for various values of 
oi < 0. The results produced by this method agree extremely well with those found with the 
direct integral equation method, as we will see in the next section. 



11 Comparing the Runge-Kutta and the direct method 

It is interesting to compare the two numerical methods used to solve the coupled set of integral 
equations. The Runge-Kutta method is a local method, which computes the function values at 
each point using the function values at neighbouring points, starting from a momentum value 
deep in the infrared region and the asymptotic expansion at that point, while the direct integral 
equation method is a global method, the complete momentum range being solved simultaneously. 
Each method employs a different set of parameters. For a given A, the Runge-Kutta method 
uses the infrared coefficients Aq and oi, while the direct method uses ^ind Fi. To compare 
results, we first have to determine the parameter sets corresponding to the same solution in the 
three-dimensional space of solutions. We run the Runge-Kutta method with X = 1, Aq = 1, and 
let ai vary till we find the solution yielding F{1) = 0.1. As mentioned before this is found for 
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ai = —10.27685. We then compute the solutions of the Runge-Kutta method at the N values 
of external momenta used in the direct integral equation method and compare the numerical 
values found with both methods, using the maximum norm. For = 81, we find 



^dir 


_pRK 


N-1 
= max 
i=0 






= 5.7 X 10"^ 




_^RK 


N-1 

= max 

i=0 






= 6.0 X 10"^ 



The agreement between these two very different numerical solution methods by far surpasses our 
initial expectations. Especially for the direct method it was hoped that the accuracy would be 
between 1/100 and 1/1000. However, the above mentioned numbers show that also this method 
achieves an even better accuracy. 

The Newton iteration of the direct method requires about 4 iterations to converge and the 
program needs approximatively 19 sec. real time to run on a Linux operated Pentium 200MHz 
PC. The Runge-Kutta method runs in approximatively 9 sec. using the Mathematica 3.0 routine 
NDSolve on the same computer. The use of two different methods is extremely important, to 
check the validity and accuracy of the solutions, especially in the case where the family of 
solutions is quite intricate. 

Although the Runge-Kutta method is faster and very accurate, it can only be used if the integral 
equations can be transformed into differential equations + boundary conditions. It also requires 
a very accurate evaluation of the starting values of the functions using the infrared asymptotic 
expansion. When the problem cannot be turned into differential equations, only the direct 
method will be usable. 

12 Including the gluon loop 

We will now briefly discuss Eqs. ^), i.e. the equations where both gluon loop and ghost loop 
are included in the gluon equation. Although it is this specific truncation which attracted our 
attention when we started the investigation of the coupled gluon- ghost equations, our physical 
expectations were better met by omitting the gluon loop. The invariances we have discussed 
yielded an unambiguous running coupling, determined by one physically relevant parameter, 
^QCD- In the following subsections, we will briefly show what changes occur when we do 
include the gluon loop and why an ambiguity occurs. 
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12.1 Symmetries of the equations 

We can repeat the analysis of Sect. ^ in the truncation we are considering now. It is easy to 
see that the solution space will still be invariant under scaling of momentum (|l3|), i.e. when 
scaling the momentum of any solution of the equations, we retrieve another solution of the same 



equations. However, the two-parameter scaling invariance (11, 12), with respect to the functions 
themselves, is now reduced to a one-parameter scaling invariance because of the additional 
constraint a = b on the scaling factors, which comes from adding the gluon loop. While the 
ghost-loop- only case was solely a function of products F{x)G{y)G{z), for various combinations 
of x,y,z, the current truncation depends on F{x)F{y)F{z) as well as on F[x)G{y)G{z). The 
fact that the three-dimensional space of solutions has lost part of its symmetry is important, 
as it means that XF{x)G'^{x) is not unique, even after an appropriate scaling of momentum. 
Globally we can say that \F{x)G'^{x) is no longer invariant, because of the admixture of \F^{x) 
terms. 

12.2 Infrared behaviour 

Because the ghost equation remains unchanged, it is easy to see that the leading infrared be- 
haviour in this case will be the same as in the ghost-loop- only case. The additional gluon loop 
in the gluon equation only yields higher order corrections. The asymptotic expansion set up in 
Sect. ^ is still generated in this case, but at some higher order it will have to be supplemented 
by other higher order series, which will be related to the leading asymptotic series. We also 
note that the power solution will not be an exact solution of the equations any more, although 
it remains the correct leading infrared asymptotic behaviour. 

12.3 Ultraviolet behaviour 

We will show that the leading log ultraviolet behaviour of the running coupling still has the 
1/logx behaviour, as expected from perturbation theory, but that the /3-coefficient is different 
from the perturbative one. This discrepancy is a bit surprising, since one expects the perturba- 
tive result to be contained in the ghost and gluon equations considered. The reason why this 
happens is that, for some reason, the pure perturbative result does not consistently solve the 
non-perturbative equations. 
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As in Sect. ^, we try the following ultraviolet solutions for F{x) and G{x), taking on the values 
Ffj^ and at some momentum /i in the perturbative regime: 



F{x) = 
G{x) = G^ 



ojlog ( - I +1 



ujIor ( - ) + 1 



n<5 



(83) 
(84) 



We check the consistency of these ultraviolet solutions by substituting the expressions in Eqs. (^, 
^). For the ghost equation Eq. (^ the treatment is identical to that of Sect. |6| and we again 
have 

7 + 25 = -l, (85) 

and 

AZiF^G^ =^(7 + 1). (86) 



Substituting the solutions Eqs. (83, 84) in the gluon equation Eq. (|^ and keeping only the 
leading log terms, we now find 



F. 



-1 







-7 


_.log 


+ 1 











F, 



-1 



-IXZiF] 



2 /" ^ 

Jx y 



wlog( -j +1 
-,27 

L^log(^) +1 



(87) 



+ XZiG 



2 z'" 

X 2y 



w log ( J ) + 1 



-,25 



After evaluation of the integrals and substitution of Eq. ([85|), 



F. 



-1 



ojlog ( - ) +1 



cjlog ( - ) +1 



7^Z,F^ 
w(27 + 1) 



wlogl -1+1 



a; log ( - ) +1 



27+1 



culog ( - I + 1 



27+1 ■ 



L^log ( - ) +1 



Consistency of this equation requires 7 < —1/3, in order to equate the leading log terms on 
both sides of the equation. We first consider the case 7 < —1/3, for which the gluon loop does 
not contribute to leading log. Then, the consistency of Eq. (Isl) requires that 



(89) 



From Eqs. (36, 89) we then find 



1 

7 = 77 



1 

3' 



(90) 
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which is inconsistent with the initial assumption 7 < —1/3. The only possibility left is 



7 = -1/3, 



(91) 



for which both the gluon and the ghost loop contribute to leading order. From Eq. ( |85D we then 
also find 

5 = -1/3, (92) 
and the condition Eq. (p6|) derived from the ghost equation yields 



27 ~ 9 
a; = -XZ^F.Gl 



(93) 



Eq. (^ then becomes 



cjlogi -j +1 

A 
+- 

UJ 



1/3 



F- 



^^log(-)+l 



nl/3 



(94) 



a 



wlog( - ) +1 



1/3 



cjlog( -) +1 



1/3 ~ 



which leads to the condition: 



a; = A 



2\Z^Fl 



3 ~ 
2 



ZxF^iG'^^ 



Eqs. (93, 95) give us 



28Zi 2 



(95) 



(96) 



llZi ^ 

which is a relation between the leading- log renormalized values of and G^, when the renor- 
malization scale [i is in the perturbative regime, in which the leading log dominates. This might 
seem to be in contradiction to perturbation theory, where the values of the renormalized quan- 
tities can take an arbitrary value and are usually fixed to 1. However, Eq. (|9^) still contains the 
renormalization constants Z\ and Z\. Taylor has shown that Z\ = \ m. the Landau gauge 
but one could still hope to be able to achieve the arbitrary renormalization of F and G by a 
suitable choice of Z\. 



If we write the far UV behaviour of F{x) and G{x) as 

F(x) ~ G log-^/^ X and G{x) ~ D log"^/^ x , 



(97) 



then, from Eqs. (|83| 84, 93, ^), the log-coefficients C and D of F{x) and Gix) are given by 

1 /7AZ. 



3Vn 



(98) 
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and 



(99) 



It is interesting to note that these leading log-coefficients are independent of the values and 
Let us now look at the ultraviolet behaviour of the running coupling. Using the solutions 



Eqs. (|83| , 34) and substituting Eq. (^), we find 



aF{x)G\x) = ^ \ . . (100) 

fAZiF.G^log ^ +1 



Now, divide numerator and denominator by XFuG^ : 



aF(x)G\x) = ,,,, — , (101) 



27 - 
4 



which can be written in the famihar form 



A'jr 

aF{x)G\x) = y ^ . (102) 

/?olog 





The leading-log coefficient Pq = 27/4 if Zi = 1, which is not in agreement with perturbation 
theory, for which (3q = 11. This seems somewhat puzzling, because all the perturbative ingredi- 
ents are contained in the non-perturbative equations, and the leading-log perturbative result can 
be retrieved from a perturbative expansion of our truncated set of Dyson-Schwinger equations. 
Nevertheless these perturbative solutions are not consistent ultraviolet asymptotic solutions of 
the non-perturbative equations themselves. The blame for this has to be put on the specific 
truncation, which loses physical information about the gauge theory. This is in contrast to the 
ghost-loop- only truncation, where the non-perturbative /5-coefficient is identical to the pertur- 
bative one, when only including ghost loops in the full gluon propagator, and which ultimately 
seems to have more physical relevance. 

12.4 Perturbative expansion 

Because of the previous remark, it is interesting to have a look at the leading order perturbative 
expansion of the truncated set of renormalized equations. To perform the expansion we set the 
zeroth order values as Fq{x) = F^ and Gq{x) = and substitute these constant values of F 
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and G in the integrals of Eq. (^) and Eq. All the integrals can be solved analytically and 
yield (for oo) 

F 

p(^x) = — (103) 

1 + A (7ZiF3 - \ZiF^Gj) log(x/M) ' 

and analogously for the ghost form factor 

Gix) = „ , ^^^^ . (104) 

l + |AZiF^G2log(x//i) 

Combining these expressions and expanding the denominator to leading order we obtain 

o dFii G?. 

aF{x)G^{x) = ^ . (105) 

1 + A (7ZiF3 + AZiF^Gl) log(x//x) 

We note that taking F^ = G^ and Zi = Zi = 1 in the previous equation would produce the 
perturbative result with /?o = 11- However, consistency with the non-perturbative integral 
equation does not allow such a choice, because of Eq. (p6|). We can even check that, using the 
condition Eq. ( p^ ) in Eq. (|105| ), we retrieve the coefficient /?o = 27/4, which means that using 
Eq. ( |96|) we find a perturbative expansion in agreement with the non-perturbative equation, 
although it is not in agreement with the usual result of pure perturbation theory. This is an 
interesting observation which deserves future research. 

12.5 Results 

We solved Eqs. (0, ^) with A = 1 and Zi = Zi = 1, for widely varying values of the parameters 
A and Fi in order to scan the two-parameter space of solutions for a given A. However, the 
loss of symmetry seems to cut out part of the solution space. Although we have made a rather 
thorough investigation of this, we will not swamp the paper with a detailed discussion since we 
think that this loss of symmetry is unphysical, thus making this truncation less interesting than 
the ghost-loop- only truncation. To put it briefiy, the fact that the ultraviolet behaviours of F{x) 
and G{x) are independent of F{ii) and G(^) destroys the invariance of the running coupling 
with respect to the choice of the individual renormalizations of F and G. Hence, choosing -F(^) 
too large will prohibit the construction of a consistent solution having the correct ultraviolet 
asymptotic behaviour. For -F(^) small this is not an obstacle, as we can show that there is an 
intermediate regime, where the log of momentum takes a different power, which allows us to 
connect to the correct ultraviolet behaviour. 
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Figure 6: Gluon form factor F{x) versus momentum x (on log- log plot), for A = 1 and Fi = 
0.001(a), 0.01(b), 0.1(c), 0.3(d) and 0.5(e). 

To see this we plot F{x) for A = 1 and Fi = 0.001, 0.01, 0.1, 0.3 and 0.5 in Fig. |. It is 
clear that the gluon form factor, which starts off as a power Ax'^'^, bends over at the cross-over 
point X, such that the further logarithmic behaviour of the function leads to a value -Fi at the 
renormalization scale x = 1. From this plot it is however clear that the curves (d, e) have a quite 
different behaviour from the others. Their ultraviolet behaviour is consistent with the log^^^^ 



analytic prediction from Sect. 12.3, while the other curves seem to show a logarithmic increase 
instead. This is of course plausible, as it is possible that the ultraviolet asymptotic behaviour 
only sets in at much higher momenta, and that in between the infrared and ultraviolet asymptotic 
behaviours there is a intermediate regime. 

From a careful investigation of the equations, we can even find a consistent analytical description 
of the intermediate regime, connecting the region of confinement to that of asymptotic freedom, 
which fits the numerical results extremely well. Consider a case where |-F(o")| <^ \G{a)\. Then, 
in the intermediate region, F{x)G'^{x) ^ F'^{x), and the gluon loop will be negligible compared 
to the ghost loop, in the gluon equation, Eq. (|7|). Keeping in mind the treatment of Sec. ^, 
we know that this has a consistent ultraviolet solution F{x) ~ log^^^ x and G{x) ~ log^^^^^ x, 
which remains valid all the way down to the region where the power behaviour bends over 
to a logarithmic behaviour. Comparison with the numerical results shows that indeed the 
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Figure 7: Running coupling a{x) versus momentum x (on log- log plot), for A = 1 and Fi 
0.001(a), 0.01(b), 0.1(c), 0.3(d) and 0.5(e). 



intermediate regime is very well reproduced by these powers of log. The ultimate ultraviolet 



behaviour of Sect. 12.3 will only set in at extremely high momentum, after the intermediate 
regime has allowed the form factors to evolve sufficiently in order to connect to the stringently 
constrained ultraviolet asymptotic behaviour. The connection of the asymptotic infrared regime, 
the intermediate log behaviour and the asymptotic ultraviolet behaviour reproduce the numerical 
result to a good accuracy. 

To evaluate the physical relevance of this truncation, it is most interesting to plot the running 
coupling in Fig. |^. We see that, in contrast to the ghost-loop- only truncation of Fig. ^, the 
various curves for the running coupling are no longer mere translations of each other on log-log 
scale; and thus, if we choose the units on each curve such that a(/i) = a^^'', we will find couplings 
which run in different ways in the intermediate regime. This means that the determination of 
the running of the strong coupling cannot be determined unambiguously in this case. 



13 Conclusions 

Following the study of von Smekal et al.0, where these authors studied the coupled system of 
Dyson-Schwinger equations for the gluon and ghost propagators, using a Ball-Chiu vertex Ansatz 
for the triple gluon vertex and a Slavnov- Taylor improved form for the gluon-gluon-ghost vertex, 
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we performed a detailed analytical and numerical analysis of the coupled gluon-ghost equations 
using the bare triple gluon and gluon-gluon-ghost vertices. The reasons to go back to the leading- 
order perturbative vertices were to avoid the ad-hoc approximations von Smekal et al. had to 
introduce in their integral equations in order to avoid logarithmic infrared singularities and to 
obtain a clear understanding of the mechanism that is the source of the new qualitative behaviour 
of the non-perturbative gluon and ghost propagators and of the running coupling. We believe 
that doing this has proven to be extremely fruitful. Firstly, the qualitative changes to the infrared 
behaviours of the propagators are solely due to the coupling of both propagator equations, and 
the details of the vertices seem to introduce merely quantitative changes. Secondly, the use of the 
bare vertices ensures that no infrared singularities occur; hence no additional approximations, 
except for the vertex Ansatze and the y-max approximation, are in principle needed in order to 
solve these equations. 

However, we did apply one more, physically motivated, truncation to the coupled gluon-ghost 
equations. Prom an analysis of the symmetries of the equations and their solutions, we believe 
that removing the gluon loop, and keeping only the ghost loop in the gluon equation, leads to a set 
of equations which is physically more relevant, because it is consistent with the renormalization 
group invar iance of the running coupling, while this is not the case in the presence of the gluon 
loop (and the approximations employed). 

We performed a detailed analytical and numerical study of the equations with and without the 
gluon loop. In the case where we removed the gluon loop, we computed the analytical asymptotic 
infrared expansion, and showed that it depends on three independent parameters defining the 
infrared behaviour of a three-dimensional family of solutions. We also derived the analytic 
ultraviolet asymptotic behaviour of the solutions, which are proportional to powers of logarithms. 
We then computed the solutions for F{x), G{x) and a{x) over the whole momentum range with 
two different numerical techniques, the Rungc-Kiitta method on the set of differential equations 
derived from the integral equations on the one hand, and on the other hand, the direct solution 
of the integral equations using a Newton iteration method to find Chebyshev approximations to 
the unknown functions. The numerical results agree very well with both asymptotic behaviours 
in the infrared and ultraviolet regions. Furthermore the results of the Runge-Kutta method 
and of the direct integral equation method agree to a very high accuracy. We found that the 
equations possess a three-dimensional family of solutions and that they all correspond to one 
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and the same physical running couphng a{x) = \F{x)G'^{x). The non-perturbative running 
couphng can be matched unambiguously to physical reality, and we showed how Aqcd can be 
determined, at least in a formal way, even though the /3-coefficient is unphysical. 

We repeated the study with inclusion of the gluon loop, and showed that this truncation is 
physically less relevant and that the non-perturbative running coupling cannot be determined 
unambiguously if one uses a bare triple gluon vertex and takes Zi to be a constant. 

To improve on the current study, we could try to reincorporate the gluon loop in the gluon 
equation in a way that respects the physical invariances of the problem. For this, we believe that 
the bare triple gluon vertex will have to be replaced by an improved vertex, like the Ball-Chiu 
vertex, and the renormalization constant Zi will have to be chosen appropriately. Furthermore 
it would be interesting to investigate the importance of the y-max approximation. 
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